3  Criterion B Details

This notebook provides step-by-step details on {term}Criterion B calculations.

4 Setup

Import Python modules

import os

import ee
import geopandas as gpd
from lonboard import Map, PolygonLayer, BitmapTileLayer

from gee_redlist.ee_auth import initialize_ee
from gee_redlist.ee_rle import load_yaml, make_eoo, area_km2, get_aoo_grid_projection

Initialize Earth Engine

initialize_ee(project=os.environ['GOOGLE_CLOUD_PROJECT'])

5 Analysis

# Define the ecosystem that this notebook is analyzing
ecosystem_code = 'MMR-T1.1.1'

# Load the country config
country_config_path = os.environ['PIXI_PROJECT_ROOT'] + '/config/country_config.yaml'
country_config = load_yaml(country_config_path)

# Extract the GEE project path from the country config
gee_project_path = country_config['gee_project_path']

# Extract the class info for the ecosystem
class_info = [x for x in country_config['classified_image']['classes'] if x['code'] == ecosystem_code][0]
print(f'{class_info = }')

ecosystem_image = {
    'asset_id': f"{gee_project_path}/{ecosystem_code}/{class_info['id']}",
    'pixel_value': class_info['id']
}
print(f'{ecosystem_image = }')

class_info = {'id': 52, 'name': 'Tanintharyi island rainforest', 'code': 'MMR-T1.1.1'}
ecosystem_image = {'asset_id': 'projects/goog-rle-assessments/assets/MMR-T1.1.1/52', 'pixel_value': 52}

classified_image_asset_id = f"{country_config['gee_project_path']}/{country_config['classified_image']['asset_id']}"
print(f'{classified_image_asset_id = }')

class_img = (
    ee.Image(classified_image_asset_id)
      .eq(ecosystem_image['pixel_value'])
      .selfMask()
)
print(f'class_img: {class_img.getInfo()}')

classified_image_asset_id = 'projects/goog-rle-assessments/assets/mm_ecosys_v7b'
class_img: {'type': 'Image', 'bands': [{'id': 'b1', 'data_type': {'type': 'PixelType', 'precision': 'int', 'min': 0, 'max': 1}, 'dimensions': [14299, 23576], 'crs': 'EPSG:4326', 'crs_transform': [0.0008084837557075694, 0, 89.60991403135986, 0, -0.0008084837557075694, 28.548369897789982]}], 'properties': {'system:footprint': {'type': 'LinearRing', 'coordinates': [[94.66768735253133, 28.548775583027798], [93.22272474215634, 28.548775580809586], [91.41652148658464, 28.548775539452144], [89.60939625454633, 28.54877430628462], [89.60949715757825, 9.485667642440385], [91.41652148658464, 9.486595817478428], [92.86148408075269, 9.486595814990817], [94.30644671824773, 9.486595784089035], [95.39016867380931, 9.48715124036664], [96.47389058147114, 9.486595770655939], [97.91885317497325, 9.486595836813008], [99.36381579356966, 9.486595799973498], [101.1708401374696, 9.48566765739743], [101.17094098108427, 28.548774348455886], [99.00257513499045, 28.548775586792935], [97.55761254104706, 28.54877552452014], [95.75140925342585, 28.548775587887942], [94.66768735253133, 28.548775583027798]]}}}

# Commented out until bug is fixed: https://github.com/developmentseed/lonboard/issues/1064
# # Determine the coordinates for viewing the image
# longitude, latitude = class_img.geometry().centroid().getInfo()['coordinates']
# longitude, latitude

6 Extent of occurrence (EOO) (subcriterion B1)

Set the scale (in meters) for reducing the image pixels to polygons. Use the image’s nominal scale unless is is less than 50 meters per pixel.

reduction_scale = max(class_img.projection().nominalScale().getInfo(), 50)
reduction_scale

90

Convert the classified image to vectors.

ecosystem_polygons = class_img.updateMask(1).reduceToVectors(
    scale=reduction_scale,
    geometry=class_img.geometry(),
    geometryType='polygon',
    maxPixels=1e12,
    bestEffort=False
)

Create a convex hull that encompasses the ecosystem polygons.

# convexHull() is called twice as a workaround for a bug
# (https://issuetracker.google.com/issues/465490917)
hull = ecosystem_polygons.geometry().convexHull(maxError=1).convexHull(maxError=1)

# Note: 
#   - area() without the projection argument calculates the geodesic area.
#   - area(proj=) calculates the planar area. the area is calculated in projected units (10 km x 10 km grid cells)!
# This is different than the javascript version, which calculates the area without reprojection
aoo_area_km2 = hull.area(1, proj=get_aoo_grid_projection()).getInfo() * 1e2
print(f'The area of the hull is {aoo_area_km2:.2f} km²')

The area of the hull is 50141.45 km²

aoo_area_km2_geodesic = hull.area().getInfo() * 1e-6
print(f'The geodesic area of the hull is {aoo_area_km2_geodesic:.2f} km²')

The geodesic area of the hull is 50336.64 km²

(aoo_area_km2_geodesic - aoo_area_km2) / aoo_area_km2

0.003892813201885356

7 TRY calculating area in projected coords

# hull_projected = hull.transform(get_aoo_grid_projection(), maxError=1)
# hull_projected.getInfo()

# hull_projected.projection().getInfo()

# hull_projected.area(1, get_aoo_grid_projection()).getInfo() * 1e2

8 DEBUG using shapely for convex hull

ecosystem_polygons_gdf = gpd.GeoDataFrame.from_features(
    ecosystem_polygons.getInfo()['features'],
    crs='EPSG:4326'
)

# Save to file
ecosystem_polygons_gdf.to_parquet('test_ecosystem_polygons.parquet')

ecosystem_polygons_gdf

	geometry	count	label
0	POLYGON ((97.70122 10.6348, 97.70203 10.6348, ...	2	1
1	POLYGON ((97.7182 12.10704, 97.71981 12.10704,...	25	1
2	POLYGON ((97.72305 12.11755, 97.72386 12.11755...	3	1
3	POLYGON ((97.72467 12.11917, 97.72547 12.11917...	275	1
4	POLYGON ((97.72952 12.13372, 97.73194 12.13372...	140	1
...	...	...	...
600	POLYGON ((97.67131 12.04722, 97.67211 12.04722...	2	1
601	POLYGON ((97.67131 12.03671, 97.67211 12.03671...	1	1
602	POLYGON ((97.67211 12.04883, 97.67373 12.04883...	3	1
603	POLYGON ((97.67373 12.03509, 97.67616 12.03509...	146	1
604	POLYGON ((97.68101 10.46259, 97.68262 10.46259...	26	1

605 rows × 3 columns

from shapely.ops import unary_union
hull_geom = unary_union(ecosystem_polygons_gdf.geometry).convex_hull

# Reproject the geometry to ESRI:54034 (World Cylindrical Equal Area) and compute area in projected units (square meters)
from pyproj import CRS, Transformer
import shapely

# Define projections
wgs84 = CRS("EPSG:4326")
cea = CRS("ESRI:54034")
project_to_cea = Transformer.from_crs(wgs84, cea, always_xy=True).transform

# Project hull geometry
hull_geom_projected = shapely.ops.transform(project_to_cea, hull_geom)

# Area is calculated in square meters (ESRI:54034's units)
hull_geom_projected.area * 1e-6

50057.86991307719

# Try to project first, then calculate hull
ecosystem_polygons_gdf_projected = ecosystem_polygons_gdf.set_crs(wgs84).to_crs(cea)
hull_geom_projected = ecosystem_polygons_gdf_projected.geometry.union_all().convex_hull
hull_geom_projected.area * 1e-6

50057.86991307719

diff = (hull_geom_projected.area  * 1e-6 - aoo_area_km2 ) / aoo_area_km2
diff

-0.0016668938119499528

9 END TEST

9.1 Display EOO Layers

9.2 Ecosystem Tiles

tile_url = class_img.getMapId(
    vis_params={
        'palette': ['blue'],
        'opacity': 0.5
    }
)['tile_fetcher'].url_format

tile_layer = BitmapTileLayer(
    data=tile_url,
    tile_size=256,
    max_requests=-1,
    min_zoom=0,
    max_zoom=19,
)

9.3 Ecosystem Polygons Layer

ecosystem_polygons_layer = PolygonLayer.from_geopandas(
    ecosystem_polygons_gdf,
    get_fill_color=[255, 0, 0, 127],
    stroked=True,
    get_line_width=2,
    get_line_color=[0, 0, 0, 150],
)

9.4 Ecosystem Hull Layer

type(hull_geom_projected)

shapely.geometry.polygon.Polygon

hull_gdf = gpd.GeoDataFrame.from_features(
    ee.FeatureCollection(hull).getInfo(),
    crs='EPSG:4326'
)

hull_layer = PolygonLayer.from_geopandas(
    hull_gdf,
    get_fill_color=[0, 0, 255, 63],
    stroked=True,
    get_line_width=200,
    get_line_color=[0, 0, 0, 150],
)

Display the map.

m = Map(
    layers=[
        ecosystem_polygons_layer,
        tile_layer,
        hull_layer,
    ],
)

m

9.5 Verify that the step-by-step results are consistent

# Direct call to `make_eoo()`
aoo_area_km2_direct_call = area_km2(make_eoo(class_img)).getInfo()
print(f'EOO area: {aoo_area_km2_direct_call:.0f} km²')

EOO area: 50337 km²

# Assert the two values are close to each other.
# assert math.isclose(aoo_area_km2, aoo_area_km2_direct_call, abs_tol=1e-4)

10 Area of Occupancy (AOO) (subcriterion B2)

The protocol for this adjustment includes the following steps:

• Intersect AOO grid with the ecosystem’s distribution map.
• Calculate extent of the ecosystem type in each grid cell (‘area’) and sum these areas to obtain the total ecosystem area (‘total area’).
• Arrange grid cells in ascending order based on their area (smaller first).
• Calculate accumulated sum of area per cell (‘cumulative area’).
• Calculate ‘cumulative proportion’ by dividing ‘cumulative area’ by ‘total area’ (cumulative proportion takes values between 0 and 1).
• Calculate AOO by counting the number of cells with a ‘cumulative proportion’ greater than 0.01 (i.e. exclude cells that in combination account for up to 1% of the total mapped extent of the ecosystem type).

10.1 Intersect AOO grid with the ecosystem’s distribution map

Load the AOO grid projection

aoo_grid_proj = get_aoo_grid_projection()

aoo_grid_proj.getInfo()

{'type': 'Projection',
 'wkt': 'PROJCS["World_Cylindrical_Equal_Area", \n  GEOGCS["WGS 84", \n    DATUM["WGS_1984", \n      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], \n      AUTHORITY["EPSG","6326"]], \n    PRIMEM["Greenwich", 0.0], \n    UNIT["degree", 0.017453292519943295], \n    AXIS["Longitude", EAST], \n    AXIS["Latitude", NORTH]], \n  PROJECTION["Cylindrical_Equal_Area"], \n  PARAMETER["central_meridian", 0.0], \n  PARAMETER["standard_parallel_1", 0.0], \n  PARAMETER["false_easting", 0.0], \n  PARAMETER["false_northing", 0.0], \n  UNIT["m", 1.0], \n  AXIS["Easting", EAST], \n  AXIS["Northing", NORTH], \n  AUTHORITY["ESRI","54034"]]',
 'transform': [10000, 0, 0, 0, 10000, 0]}

Extract the grid scale parameters

aoo_x_scale, _, _, _, aoo_y_scale, _ = aoo_grid_proj.getInfo()['transform']
print(f'{aoo_x_scale = } meters')
print(f'{aoo_y_scale = } meters')

aoo_x_scale = 10000 meters
aoo_y_scale = 10000 meters

• Create an Earth Engine feature collection of AOO grid cells that intersect with the ecosystem, and calculate the fractional coverage of the ecosystem within the grid cell.

fractional_coverage_fc = class_img.unmask().reduceRegions(
  collection=class_img.geometry().coveringGrid(aoo_grid_proj),
  reducer=ee.Reducer.mean(),
).filter(ee.Filter.gt('mean', 0))

# Convert the Earth Engine feature collection to a GeoPandas GeoDataFrame.
fractional_coverage_gdf = ee.data.computeFeatures({
    "expression": fractional_coverage_fc,
    "fileFormat": "GEOPANDAS_GEODATAFRAME",
})
fractional_coverage_gdf.rename(columns={"mean": "coverage"}, inplace=True)

# Set the CRS (Earth Engine data is in EPSG:4326)
fractional_coverage_gdf = fractional_coverage_gdf.set_crs('EPSG:4326')

aoo_grid_cell_count = len(fractional_coverage_gdf)
aoo_grid_cell_count

206

10.2 Calculate grid cell area and the total ecosystem area

• Calculate extent of the ecosystem type in each grid cell (‘area’) and sum these areas to obtain the total ecosystem area (‘total area’).

fractional_coverage_gdf['area'] = fractional_coverage_gdf['coverage'] * aoo_x_scale * aoo_y_scale

fractional_coverage_gdf.sort_values(by="area")[0:4]

	geometry	coverage	area
175	POLYGON ((98.09603 12.57845, 98.18586 12.57845...	0.000035	3481.367444
125	POLYGON ((98.09603 11.93102, 98.18586 11.93102...	0.000079	7875.324046
111	POLYGON ((98.09603 11.74632, 98.18586 11.74632...	0.000079	7880.464175
74	POLYGON ((98.54519 11.1008, 98.63502 11.1008, ...	0.000079	7898.091697

total_area_km2 = fractional_coverage_gdf['area'].sum() / 1e6
print(f'Total ecosystem area: {total_area_km2:.0f} km²')

Total ecosystem area: 1937 km²

10.3 Calculate cumulative area in ordered cells

• Arrange grid cells in ascending order based on their area (smaller first).
• Calculate accumulated sum of area per cell (‘cumulative area’).

fractional_coverage_gdf = fractional_coverage_gdf.sort_values(by="area")
fractional_coverage_gdf["cumulative_area"] = fractional_coverage_gdf["area"].cumsum()

fractional_coverage_gdf.sort_values(by="area").head()

	geometry	coverage	area	cumulative_area
175	POLYGON ((98.09603 12.57845, 98.18586 12.57845...	0.000035	3481.367444	3481.367444
125	POLYGON ((98.09603 11.93102, 98.18586 11.93102...	0.000079	7875.324046	11356.691491
111	POLYGON ((98.09603 11.74632, 98.18586 11.74632...	0.000079	7880.464175	19237.155665
74	POLYGON ((98.54519 11.1008, 98.63502 11.1008, ...	0.000079	7898.091697	27135.247363
160	POLYGON ((97.91637 12.39331, 98.0062 12.39331,...	0.000123	12301.457861	39436.705224

10.4 Calculate the cumulative proportion

• Calculate ‘cumulative proportion’ by dividing ‘cumulative area’ by ‘total area’ (cumulative proportion takes values between 0 and 1).

fractional_coverage_gdf["cumulative_proportion"] = fractional_coverage_gdf["cumulative_area"] / 1e6 / total_area_km2

fractional_coverage_gdf

	geometry	coverage	area	cumulative_area	cumulative_proportion
175	POLYGON ((98.09603 12.57845, 98.18586 12.57845...	0.000035	3.481367e+03	3.481367e+03	0.000002
125	POLYGON ((98.09603 11.93102, 98.18586 11.93102...	0.000079	7.875324e+03	1.135669e+04	0.000006
111	POLYGON ((98.09603 11.74632, 98.18586 11.74632...	0.000079	7.880464e+03	1.923716e+04	0.000010
74	POLYGON ((98.54519 11.1008, 98.63502 11.1008, ...	0.000079	7.898092e+03	2.713525e+04	0.000014
160	POLYGON ((97.91637 12.39331, 98.0062 12.39331,...	0.000123	1.230146e+04	3.943671e+04	0.000020
...	...	...	...	...	...
108	POLYGON ((98.45536 11.65401, 98.54519 11.65401...	0.549146	5.491457e+07	1.689946e+09	0.872610
102	POLYGON ((98.45536 11.56174, 98.54519 11.56174...	0.561263	5.612625e+07	1.746072e+09	0.901591
164	POLYGON ((98.36552 12.39331, 98.45536 12.39331...	0.586257	5.862569e+07	1.804698e+09	0.931863
99	POLYGON ((98.18586 11.56174, 98.27569 11.56174...	0.607082	6.070816e+07	1.865406e+09	0.963210
93	POLYGON ((98.18586 11.46949, 98.27569 11.46949...	0.712497	7.124972e+07	1.936656e+09	1.000000

206 rows × 5 columns

10.5 Calculate AOO

• Calculate AOO by counting the number of cells with a ‘cumulative proportion’ greater than 0.01 (i.e. exclude cells that in combination account for up to 1% of the total mapped extent of the ecosystem type).

aoo_grid_cells = fractional_coverage_gdf[fractional_coverage_gdf["cumulative_proportion"] > 0.01]
print(f"AOO (number of cells with cumulative proportion > 0.01): {len(aoo_grid_cells)}")

AOO (number of cells with cumulative proportion > 0.01): 139

aoo_grid_cells_dropped = fractional_coverage_gdf[fractional_coverage_gdf["cumulative_proportion"] <= 0.01]

10.6 Display the layers

aoo_grid_cells_layer = PolygonLayer.from_geopandas(
    aoo_grid_cells,
    get_fill_color=[0, 255, 0, 63],
)

aoo_grid_cells_layer_dropped = PolygonLayer.from_geopandas(
    aoo_grid_cells_dropped,
    get_fill_color=[255, 0, 0, 63],
)

m = Map(
    layers=[
        tile_layer,
        aoo_grid_cells_layer,
        aoo_grid_cells_layer_dropped
    ],
)

m

10.7 Verify that the step-by-step results are consistent

# # Direct call to `make_aoo()`
# aoo_cells_direct_call = make_aoo(class_img)

# print(f'AOO: {aoo_grid_cell_count_direct_call} grid cells')

# assert aoo_grid_cell_count == aoo_grid_cell_count_direct_call

11 Criterion B Summary

from IPython.display import Markdown, display
display(
    Markdown(f'AOO and EOO were measured as '
      f'{aoo_grid_cell_count} 10 x 10 km grid cells '
      f'and {aoo_area_km2:.0f} km², respectively. '
      f'See [Criterion B Details](mmr-t1-1-1-crit-b-test) for more information. '
      f'There is no evidence that suggests this ecosystem should meet the criteria to '
      f'be listed as Near Threatened. The ecosystem is assessed as Least Concern under '
      f'Criterion B1 and B2. '
      f'**Least Concern.**'
    )
)

AOO and EOO were measured as 206 10 x 10 km grid cells and 50141 km², respectively. See Criterion B Details for more information. There is no evidence that suggests this ecosystem should meet the criteria to be listed as Near Threatened. The ecosystem is assessed as Least Concern under Criterion B1 and B2. Least Concern.